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Abstract — Stable periodic behavior is the desired behavior 
of many robotic systems, such as walking robots, swimming 
robots, and other biologically inspired systems. Designing and 
verifying controllers for stable periodic behaviour is signifi- 
cantly more challenging than for equilibria. We study pointwise 
contraction conditions for an autonomous system which imply 
the existence and global stability of a limit cycle in a compact 
manifold. Unlike methods based on Lyapunov theory, our 
condition does not require knowledge of the location of the 
limit cycle in state-space. This allows straightforward extension 
to robustness via parameter-dependent certificates, and compo- 
sitions of systems in feedback and serial interconnections. This 
new criterion also has interesting connections to contraction 
theory for trajectories, orbital stability analysis via transverse 
linearization, and topological facts about regions of attraction 
of limit cycles. 

I. Introduction 

Systems with periodic solutions are central to biologically- 
inspired robot locomotion [1], [2], and are also in many other 
areas of science and engineering, e.g. phase-locked loops or 
vortex shedding from aircraft. In biology, oscillating systems 
seem to be the rule rather than the exception [3]. 

The basic question we address in this paper is the follow- 
ing; when does a system of the form 

i^fix) (1) 

have the property that all solutions starting from a particular 
set K converge asymptotically to a unique limit cycle? This 
is the property of orbital stability, and is substantially more 
difficult to deal with than simple stability of an equilibrium. 

There is a long and distinguished history of research into 
limit cycles for nonlinear systems. For example, the famous 
result of Poincare-Bendixson gives a very simple condition 
for planar systems. However this result relies heavily on the 
Jordan curve theorem and for higher-dimensional systems 
fails to hold: bounded orbits of smooth systems may have 
chaotic behaviour An important generalization to monotone 
cyclic feedback systems was published in [4], however this 
depends on quite a special system structure and there are 
many application areas where it does not apply. Birkhoff 
gave necessary conditions for periodic solutions in terms of 
the existence of particular "phase variables", or associated 
differential one-forms [5], [6]. 

However, all of these conditions imply the existence of at 
least one limit cycle, but give no insight into the number 



of limit cycles, or their stability. Previous work by the first 
author and colleagues has developed the Lyapunov-based 
"transverse dynamics" approach to study regions of stability 
of limit cycles [7], [8], [9], [10], however this approach is not 
applicable when the system dynamics are uncertain, since 
uncertainty will generally change the location of the limit 
cycle in state space. 

In recent years many efficient computational methods for 
proving stability of equilibria of nonlinear systems have 
been proposed, based on searching for "stability certificates" 
such as Lyapunov functions and contraction metrics. For 
the purposes of robust stability analysis, a central difference 
between these two approaches is that a Lyapunov function 
must generally be constructed about a known equilibrium, 
whereas a contraction metric implies the existence of an 
equilibrium indirectly [11]. This can be useful for robustness 
analysis of highly nonlinear systems where the equilibrium 
point may change depending on the unknown dynamics [12]. 

In this paper we extend contraction analysis to the study 
of limit cycles. It will be shown that the existence of 
a contraction metric of a particular "transverse" form is 
equivalent to the existence of a unique orbitally stable 
limit cycle. In other words, we decouple the property of 
convergence from the task of finding a particular limiting 
solution. A major benefit of this approach is that conditions 
for stability can be extended to conditions for robustness, i.e. 
existence of a stable limit cycle under parameter uncertainty 
and interconnection with other (possibly uncertain) systems. 

II. Problem Setup and Preliminaries 

We assume that / : if — >• K" in ([T]i is smooth and a; e K". 
In the situations we will consider in this paper the solution 
of ([TJ will always exist and be unique. 

A periodic solution x* is one for which there exists some 
T > such that x*{t) = + T) for all t. Equilibria ai-e 
trivially periodic for every T, but for oscillatory solutions - 
which are our main concern - there is some minimal time T 
such that the above holds and this is referred to as the period. 
The orbit of a periodic solution is the set X* :— {x : x = 
x*{t) for some t}. Note that while periodic solutions cannot 
be asymptotically stable, their orbits can be, and in this case 
we may say that the solution is orbitally stable (see, e.g.. 




< 1 



Fig. 1. Regions of attraction to limit cycles are homeomorphic to discs 
parametrized by a scalar phase variable r. 



[13]). We will also refer to this as "convergence to a periodic 
orbit". 

The basic statement is: give pointwise conditions on the 
vector field /(•) and a compact set K C M" guaranteeing 
that every solution of ([T]i starting from an initial condition 
in K converges to a unique periodic orbit X* as i oo. 

III. Stability and Regions of Attraction of 
Periodic Solutions 

In this section we discuss some aspects of oscillating 
systems that render the study of stability substantially more 
challenging that stationary systems. 

It is well known that periodic solutions of an autonomous 
differential equation can never be locally asymptotically 
stable. This is clear from the fact that solutions which have 
initial conditions on the periodic orbit but offset in time will 
never approach each other Orbital stability is the strongest 
form of stability that can be hoped for. 

There are also interesting properties of the "global" struc- 
ture of regions of attraction to periodic orbits. It is known that 
the the region of attraction is a continuous deformation of the 
set D"^^ X S^, i.e. the cartesian product of an open unit disc 
of dimension one less than the full state space, with a single 
circle coordinate. These are often referred to as "transversal" 
and "phase" coordinates, respectively. In all cases except 
possibly with n = 5 it is guaranteed that the deformation 
is smooth differentiable, due to the recent resolution of the 
Poincare Hypothesis [6]. 

In that case, let us consider such a change of variables: 



x± =p{x), 



q{x) 



(2) 



with p{x) E M"^^ and t e which we assume to be 
differentiable. For r to be considered a "phase variable" it 
must satisfy t = ^f{x) > for all x G K. 

Indeed, the existence of a global phase variable - a scalar 
variable that is always increasing - implies the existence of 
at least one periodic orbit, though not stability [5]. 

A change of variables into phase and transversal coor- 
dinates depends on the existence of a codimension-1 fo- 
liation, i.e. a family of n — 1-dimensional submanifolds 



("leaves") which are parametrized by the phase variable. 
Not all nonsingular vector fields admit such a foliation, e.g. 
Hopf bifurcations do not [14]. Regardless, even if such a 
transformation exists it may be prohibitively hard to find. 

Nevertheless, let us mention briefly some important points. 
From the definition of the transformations, we have 

dp 



and 



dx 
dq 



{x)f{x), 



dx 



{x)f{x), 



noting that f is a scalar greater than zero. 

It is clear that r can be considered a reparametrization of 
time, and then one can analyze instead of ij^. Since 
T e §^ the differential equation 

dx±^ 1 dp ^ 



-{x)f{x) =: /j_(a;_L,r) 



can be considered to be a periodically time-varying differ- 
ential equation of dimension n — 1. Now, it is known that 
systems which are periodically time-varying and contracting 
have the property that all solutions converge to a particular 
periodic solution [11]. Making the "lifting" transformation 
back to considering r to be a state, this then implies 
convergence to a periodic orbit in dimension n. Thus if one 
can construct a transformation to phase and transversal coor- 
dinates, and prove that the transversal system is contracting, 
then one has proven orbital stability. 

The remarkable fact is that one can show transverse 
contraction without ever constructing the global change of 
variables. 

IV. Contraction Conditions for Limit Cycles 

The following theorem relates existence of a special form 
of metric to a property sometimes referred to as Zhukovsky 
stability, which is deeply linked to the existence of periodic 
solutions [15]. 

Theorem 1: Let K C M" be compact, path-connected, 
and strictly forward invariant. If there exists a smooth 
positive-definite matrix function M{x) on K and a A > 
such that 



S'{M + MF + F'M + \M)5 < 0, 



(3) 



for all S e M" satisfying S'M{x)f{x) = 0, then for 
every two solutions xi and X2 with initial conditions in K 
there exists time reparametrizations Ti{t), T2{t) such that 

Xi{Ti{t)) -> X2[T2{t)) as t ^ CO. 

Proof: since the set K is path-connected by definition 
there exists a path between any two points xq and xi in 
K. Such a path can be considered as a smooth mapping 
7 : [0, 1] -> M" with 7(0) = xq and 7(1) = xi. 

Denote by r(a;o,a;i) the set of all such smooth paths 
between xq and xi remaining in K and associate with each 
a length 



d 



7(s) 



M(7(s))^7(s)ds 



and then the standard Riemannian metric between two points 
xq and is 

d{xo,xi)= inf L(7) 

76r(a;o,2;i) 

As in contraction analysis, we study the instantaneous 
change in the length ^(7) as xi and X2 evolve along 
solutions of i; = f{x): 

|(|;7(.))'m(7W)|7(.) 

However, since we are allowing time reparametrizations 
between solutions, we also allow each point on 7(5) to move 
at a different rate ^^-^r^ = /(7(s))a(7(s)). Then 

Without loss of generality for this pointwise condition we 
can let a{s) = 1 and denote z :— ^^^^^j^, then in terms of 
the virtual displacement 6 we have 



d 
dt 



S = F{x)5 + zf{x) 

{5' MS) =(5' {f{x)'M{x) + M{x)F{x) + A/(x)) 5 
+ 2z5'M{x)f{x) 



Now, since z is a free variable, there always exists a z mak- 
ing the evolution of this length segment negatative, except 
possibly when 5'M{x)f{x) = and z has no influence. 
From this we derive the transverse contraction condition in 
the theorem. □ 

It is known that systems satisfying such a property have 
limit cycles [15], but with the framework of contraction the 
proofs are simpler, so we give a proof here. 

Theorem 2: If the conditions of Theorem 1 are satisfied, 
then all solutions starting with x(0) e K converge to a 
unique limit cycle. 

Proof: since K is invariant and compact, it follows that 
0,{xo) exists and is a compact subset of K. Furthermore, a 
clear implication of Theorem [T] is that all points in K have 
the same cj-limit set, which we denote Vl{K). 

Pick a point x* in Q,{K). By strict invariance, this is 
an interior point of K. By assumption, f{x*) 7^ so we 
can construct the hyperplane orthogonal to f{x*), which we 
denote S. We will prove convergence to a limit cycle by 
constructing a Poincare map on S. 

Since /(•) is smooth, for x in some neighborhood B of x* 
we have that f{xyf{x*) > 0, so in Bs '.^ B n S solution 
curves are transversal to S and pass through it in the same 
direction as at x*. 

Since x* is in the cj-limit set for all points in K, and Bg 
is transversal, the evolution of the system from any point 
xo{t) S Bs eventually passes through ^5 again, i.e. xo{t + 
s) S Bs where s > depends on xq. This evolution can be 
represented by a Poincare map T : Bs — > Bs- 

Take the distance between two points d{xo,xi) on Bs 
to be Riemannian metric distance from Lemma 1 . Note that 



although the two points lie on the n — 1 dimensional set Bs, 
the curves joining them in the definition of d may pass out 
of the plane and through n-dimensional space. By Lemma 
1, we have that d{T{xo),T{xi)) < d{xo,xi). Hence T is 
a contractive map from Bs onto itself, and by the Banach 
fixed point theorem has a unique stable fixed point, which 
is its only limit point so must be x*. By standard results on 
Poincare maps this implies that x* is a point on a limit cycle, 
which by Proposition 1 all solutions converge to. 
□ 

The new conditions is weaker than contraction, which 
implies that any two solutions converge, i.e. xo{t) xi{t), 
but stronger than orbital stability, which implies only conver- 
gence to the same w-limit set and thus doesn't preclude chaos 
[15]. Our condition implies a behaviour that is sometimes 
referred to as Zhukovsky stability and has been used in 
several recent papers on limit cycle stability and system 
identification, see e.g. [16], [9], [17] and apparently goes 
back to Poincare in its essential argument [13]. 

V. Singular Contraction Metrics 

The following equivalent condition is easier to use in 
some applications, where a set of transversal coordinates are 
available [9]. 

Theorem 3: Given a particular choice of normal vector 
z{x) satisfying z{x)' f{x) > every wheraM and a particular 
set of transversal coordinates, condition Q is equivalent to 
the existence of a singular positive-semidefinite contraction 
metric Ms{x) with rank[Afs(a;)] = n—1 and transversal, i.e. 
Ms{x)z{x) = 0, which satisfies 

Ms + MsF + F'Ms + AM, < 0. (4) 
Proof: Let tt{x) G M^^" be the orthogonal projection 
matrix on to the subspace spanned by z{x). Let n(a;) G 
jj(n-i)xn jjjg projection on to the subspace orthogonal 
to z{x), i.e. on to the transversal coordinates. Let f± = 
Il{x)f{x) and /y = Tr{x)f(x) be the projections of f{x) on 
to these subspaces. 

Decompose the singular metric like so: 

Ms{x) = n(x)'M(a;)n(a;), 

then we can construct a new positive-definite metric M{x) 
in the following way 



M(x) 



-U{x) 


1 r 


t:{x) 





M{x) 



1 





"n(x)" 




t:{x) 



The interpretation of this is that a differential 5 is decom- 
posed into a component collinear with f{x) which we can 
denote by (5|[ and a component orthogonal to this denoted by 
5^. Then 

5'M{x)5 = 5'^M{x)5^ - 26'^M f^5\\/ f\\ + \5\\\'^. 

Note that f\\ ^ since j(x) ^ and j{x)'z{x) > 0. 
Now, the conditions of Theorem[T|are that S'M{x)f{x) = 0. 

'Taking z{x) = f{x)/\f{x)\ is one possible choice, but in some 
applications freedom in z{x) is beneficial. 



By the above construction of M{x), the terms involving 6± 
cancel and so S'M{x)f(x) — implies 

5||(-/iM/i + /|f) =0 

But since Mg is positive definite, and f± and /|| are projec- 
tions on to an orthonormal basis of f{x), and f{x) 7^ 0, it 
follows that the term in the brackets can not be zero, so this 
implies that = 0. This implies that 

±i6'Mix)S) = j^i6'Mx)S^) 

which by Condition (|4]) is negative, and hence the system 
satisfies the conditions of Theorem [T] 

Given a metric M{x) satisfying the conditions of The- 
orem [T| there exists a stable periodic solution which all 
initial conditions converge to. By constructing a differential 
system explicitly in the transversal coordinates defined be 
the condition z{xyS = using the method described in [9] 
one has a stable time-varying linear system of dimension 
n — 1. By standard results on Lyapunov functions for time- 
varying Unear systems, along each time-reversed trajectory 
there exists a smooth solution to the Lyapunov differential 
equation, which can be used as a singular contraction metric. 
See [9] for more details on such constructions. 
□ 

This formulation is also easy to extend to hybrid dynam- 
ical systems (e.g. walking robots) by enforcing a discrete 
contraction condition on switching surfaces, see [11] and[9] 
for details. 

A. Generalized Jacobian and Transverse Linearization 

In this section we show given the existence of a singular 
contraction metric, the negative definiteness of an associated 
generalized Jacobian is equivalent to contraction of transver- 
sal coordinates. 

Let 

M,{x) = n{xyM{x)u{x) 

where M{x) E M"^^ and is positive definite. n(a;) : M" — > 
R""^ is a projection on to the subspace orthogonal to f{x). 

Consider a decomposition of M{x) — L{x)' L{x) with 
L{x) nonsingular, then set Q{x) := L{x)Il{x) and consider 
the generalized Jacobian of the system under the change of 
coordinates z — L{x)Il{x)x 

where ji denotes pseudoinverse. Note that 9" = 
Il{xy L{x)~^ . Then we have 

F = Lun'L-^ + L{tiu' + nFn')L-\ 

The term nil' + IIFII' is exactly the transverse linearization 
from [9] with orthogonal hypersurfaces, which we denote 
F±. In fact, in the trivial case in which the identity is a 
contraction metric, the generalized Jacobian is the transverse 
linearization. Now, since 11 is a projection on to a lower 
dimensional space, IIII' = /, so the above can be simplified 
as 

F = LL-^ +LFj_L-^ 



If the generalized Jacobian is uniformly negative: 

F + F < -\I 

for some A > 0, we can multiply both sides of the above 
inequality by L and obtain 

L'L + L'L + L'LF±_ + F'^L' L < -XL' L 

And since M ~ L'L this is 

M + MFj_ +F'j^M < \M 

which is exactly the conditions for contraction of the trans- 
verse dynamics [11]. Furthermore, when a periodic solution 
is known, then this is also a Lyapunov ODE condition, which 
is satisfied if and only if the periodic solution is orbitally 
exponentially stable [9]. 

Note that in [9] the transverse dynamics were computed 
explicitly about a known periodic solution and used to 
compute regions of orbital stability. In this paper, we consider 
the linearization of the transverse dynamics about every 
solution of the system, though only implicitly. In general, 
it is difficult to compute an analytical formula for L(a;): 
in principle it could be done via the Cholesky algorithm 
although the expressions would be extremely complicated. 

Given the results of [11], the following theorem is an easy 
consequence of the above derivation: 

Theorem 4: If there exists a 11 (x), L so that F is con- 
tracting with respect to some matrix measure on a compact 
forward-invariant set, then there exists a unique orbitally 
stable limit cycle. 

VI. Extensions of the Basic Theory 

The results in the previous sections are derived for au- 
tonomous systems, i.e. time-invariant systems with no exoge- 
nous inputs. However, the transverse contraction machinery 
can be applied for a much wider range of systems including 
compositions, feedback interconnections, and systems with 
bounded disturbances. 

In many applications in which exact model parameters 
are unavailable, it is desirable to characterize parameter 
ranges over which the qualitative behaviour of the system 
remains the same. Engineering motivations are well known, 
but robustness analysis has also become of interest recently 
in biology, including as a measure of model validity [23]. In 
particular, the papers [23], [24] address robustness of limit 
cycles by gridding over parameter ranges and simulating the 
nonlinear system until convergence can be ascertained. Grid- 
ding and simulation becomes very expensive computationally 
as the number of parameters grows, so alternative methods 
are desirable. 

Robustness analysis of equilibria and trajectories is clas- 
sical topic for which many methods have been derived, but 
robustness of limit cycles is substantially more challenging 
to address. In [25] small-gain conditions were derived for 
limit-cycle robustness, though such conditions are often very 
conservative. 



A. Feedback Interconnections and Incremental Dissipativity 
For this result we consider systems with exogenous inputs: 

x = f{x,w) (5) 

We define a system as transverse incrementally dissipative 
(TID) with a supply rate a{x^'w,5x,5w) if the following 
condition holds 

X 7 ) ■ 

for all dx such that S'^M{x)f{x,w) = 0. Note that if 
w — then the above condition reduces to non-strict 
transverse contraction. Hence, every TID system is transverse 
contracting. 

If the system Q also has an output y = g{x,w) with 
Jacobian ^ = G{x,'w), then we can mention two partic- 
ularly useful informative supply rates: incremental versions 
of classical small-gain: 



0-7(4, = - \G{x,w)Sx\'' 



and passivity: 



CTpiSx.Sw) = S'^G{x,w)Sa;, 



where passivity requires that input w and output y have the 
same dimension. 

We now give the following theorem: 

Theorems.- Given two systems xi — fi{xi,wi) and 
2^2 — f2{x2,W2), which are interconnected so that wi — 
52(3^2) and W2 = gi{xi). Then if the two systems are TID 
with supply rates cri,a2 such that ai + a2 < — a|(52:P for 
some a > 0, then the interconnected system is transverse 
contracting. 

Proof Let Si and 62 refer to the state differentials of system 
1 and system 2, respectively, and — [6162]' ■ Furthermore, 
let f{x) = [fi{xi,g2{x2)) .f2{x2,9i{xi))]' be the dynamics 
of the interconnected system. 

Suppose the two systems are verified to be TID with 
metrics Afi(a;i) and M2{x2), then consider the metric for 
the interconnected system 



M = 



Ah{xi) 






M2(X2) 



It is clear from construction that this satisfies 
d 



dt 



for all Sx satisfying 5'^M{x)f{x) — and for some A > 0, 
where the second inequality follows from the compactness 
of K and continuity and positive-definiteness of M{x). □. 

The above method is related to the so-called S-Procedure 
in robust control, widely used in H°° and Integral Quadratic 
Constraint analysis [18], [19]. 

Note that TID is a stronger property than incremental 
dissipativity so, for example, any linear passive system 
satisfies the TID property with supply rate dp, and any linear 
stable system satisfies TID with supply rate where 7 is 
the square of the H°° norm of the system (peak frequency 



response). The interconnections of oscillators with passive 
systems has been studied in the past via describing function 
methods [20], but without rigorous proofs of stability. 

This formulation may also be helpful in analyzing inter- 
connections of oscillators [21], [22]. Note that we have a 
condition by which two non-identical oscillators will, when 
interconnected, converge to a periodic solution which may 
be entirely different than their natural oscillations. 

B. Serial Compositions of Systems 

A relatively simple application of the above theorem is 
to consider the composition of a contracting system and a 
transverse-contracting system. 



Xi 
X2 



,f2ixi,X2) 



(6) 
(7) 



Theorem 6: Suppose for each fixed xi, /2 is transverse 
contracting with metric Al2{x2,xi), i.e. 

S'2iF2M22 + M2F22 + 7^M2/(X2) + A2M2)4 < 
0X2 

for all 62 satisfying 62^12/2 = and /i is "fully" contracting 
in the sense of [11], i.e. there exists Mi{xi) such that 

F[Mi + MiFi + ^MJ{xi) + AiMi < 
ax I 

then the composed system is transverse contracting, and 
hence has a unique stable limit cycle. 

Proof we prove this theorem by constructing a metric which 
decomposes as 

6'Mc5 5[Mi5i + aS'2M2S2 = 

which will be shown to verify the existence of a unique 
stable limit cycle. Let x = \x\x'^ and /(.t) = 
[/i(a;i)'/2(a;i, ^2)']'- Since the contraction conditions are 
homogeneous with respect to b, and (5 = is trivial, it is 
sufficient to consider the case where \5\ = 1. 

First, we note that since /2(a;) 7^ is X and K is compact, 
there exists an e > such that |/2(2:)| > e for all x K. 

Second, since system 1 is contracting and K is compact, 
/i uniformly [11]. The transversality condition for the 
metric Mc is 

5[Mih+5'2M2f2=0 (8) 

However, since /2 is bounded below and fi converges 
uniformly to zero, the normal vector to the surface defined 
by ([8]) converges to that defined by ^2^2/2 = and hence 
the compact sets of 5 of norm one satisfying these conditions 
converge uniformly. 

Now, d/ dt[5' McS\ = 6'HS where H decomposes into 
blocks corresponding to 61 and 62 like so: 



H = 



a{F{Mi + MiFi + Ml) 
F^iAfi +Mii^2i 



F^^M2 + MiF2i 
F^^M2 + M2F22 + M2 

Consider the fixed x\ to which the contracting system xi = 
fi{xi) converges, so that fi{x* — 0). Now, since system 1 



is contracting, the upper-left block is negative definite and 
then by the Schur complement it follows that the maximum 
value of 5'H6 on the subspace satisfying S'Mf = can be 
made strictly negative by choosing a sufficiently large. 

Due to continuity of H and the convergence of the sets 
of 6, this implies that from any initial conditions there exists 
a finite time after which S'H6 < for all S satisfying ([HJ, 
which implies the existence of a unique stable periodic orbit 
by Theorems [T] and |2] 
□. 

The opposite composition, a transverse-contracting system 
driving a contracting system clearly converges to a periodic 
solution due to natural input-to-state stability properties of 
contracting systems. In a sense, the second system can be 
considered as being driven by a periodic input. 

C. Robustness to Parametric Variation 

Suppose the system dynamics depend on some parameter 
vector 9, i.e. 

When studying robustness of equilibria of such systems, 
a widely-used method is to search for either a parameter- 
dependent Lyapunov function via an optimization procedure, 
see, e.g., [26], [27]. 

In the context of the present paper, we assume that a 
particular set K is robustly forward invariant - which can be 
verified using the methods of [26] - then robust existence of a 
single globally stable (within K) Hmit cycle is ensured if one 
can find a parameter-dependent contraction metric M {x, 9) 
which satisfies 

5'{M{x, 9) + 2F{x, 9)'M{x, 9) + \M{x, 9))5 < 

for all S such that S'M{x,9)f{x,9) = 0, and for all xeK 
and 9 in some set 9, where i^T is a forward invariant set. 

D. Robustness to Bounded Disturbance 

Consider the global coordinates x±,t - either implicitly 
or explicitly defined. Since t E the dynamics of x± 
can be considered a periodic differential equation with a 
transformation of time. This makes it clear that any internal 
perturbation in / which keeps f > and F± (x) contracting 
still results in a limit cycle (c.f. above). 

Bounded external perturbations will also have bounded 
effect on behavior. Letting R — J^^^' S'MS we have 

R + 2XR < 

Consider a bounded external disturbance, i.e. x = f{x) + 
d{t), where \d\ < c?max, then we have 

R + 2XR<\nix)d{t)\Mi^) 

so after exponentially-forgotten transients, the perturbed sys- 
tem is within a ball of radius R. For further details on such 
analysis, see [11]. 
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